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Abstract 

A consolidated mathematical formulation of the spherically sym- 
metric mass-transfer problem is presented, with the quasi-stationary 
approximating equations derived from a perturbation point of view for 
the leading-order effect. For the diffusion-controlled quasi-stationary 
process, a mathematically complete set of the exact analytical solu- 
tions is obtained in implicit forms to cover the entire parameter range. 
Furthermore, accurate explicit formulas for the particle radius as a 
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function of time are also constructed semi-empirically for convenience 
in engineering practice. Both dissolution of a particle in a solvent 
and growth of it by precipitation in a supersaturated environment are 
considered in the present work. 
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1 Introduction 



Many technological applications involve mass transfer with respect to parti- 
cles. For fundamental understanding in mathematical terms, the problem of 
mass transfer to and from a particle is typically treated as an isolated sphere 
with time-dependent radius in a continuous medium of infinite extent as a 



consequence of heat-mass transfer (cf. Scriven, 1959 Duda and Vrentas 



1969 Cable and Frade, 1987; Rice and Do 2006). When the attention is 



focused on the dissolution (or growth by precipitation) of solid particles in 
liquids, as especially important in pharmaceutical dosage form development, 
the mass transfer problem may often be simplified by ignoring the effects of 
convection and phase-change heating such that the governing equations be- 
come linear with the "quasi-stationary" treatment. Thus, the mathematical 
problem is tractable for deriving analytical solutions as usually desired for 



engineering evaluations. Moreover, when the mass transfer is mainly limited 
by the diffusion process rather than the rate of phase change, as often to 
be the case in several realistic applications, the solute concentration at the 
particle-medium interface can be assumed to take a constant value of the 
so-called solubility. Then the mathematical problem physically describes a 
diffusion-controlled mass transfer process, with all the boundary conditions 
given in the form of Dirichlet type. Despite the efforts of many authors 
over years, the mathematical analyses of this relatively simplified diffusion- 
controlled quasi-stationary mass transfer problem have not been thoroughly 
satisfactory in terms of completeness and clarity. Basic understanding of the 
accuracy and validity of some approximation formulas seems to be lacking 
in the literature. 

The purpose here is to first present a consolidated mathematical formu- 
lation of the spherically symmetric mass-transfer problem, then to derive the 
quasi-stationary approximating equations mainly based on a perturbation 
procedure for the leading-order effect, and to provide a complete set of ex- 
act analytical solutions for the entire parameter range. Because the exact 
solutions can only be written in implicit forms, effort in semi-empirical con- 
struction of explicit formulas for the particle radius as a function of time is 
also made for convenience in engineering practice. 
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2 Problem Formulation 



The diffusion-controlled mass transfer to and from a spherical solid particle 
of (a time dependent) radius R in an incompressible continuous fluid medium 
with a constant density pm and a constant diffusion coefficient D is governed 



by (Bird et al. 1960 p. 557) 



dC dC D d l^odC' 



di 



dr dr \ dr 



R < r < oo 



(1) 



where C denotes the mass concentration of the solute (namely the dissolved 
solid from the particle), t the time, and r the radial distance from the center 
of the sphere. In an incompressible fluid with a spherically symmetric flow, 
the radial velocity is simply 



Pm, 



R' 

r . 



dR 

di 



(2) 



where pp denotes the (constant) solid particle density, to satisfy the equation 
of continuity and to account for the effect of volume change during the solute 



phase change (e.g., Scriven, 1959). At the particle surface, the mass balance 



based on Pick's first law of binary diffusion accounting for the bulk flow effect 



with the solvent flux being ignored (Bird et al. , 1960, p. 502) leads to 



dR 

di 



D 



dC 



Pp{l - C/pra) dr 



(3) 



f=R 
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In a diffusion-controlled process, the typical boundary conditions for C 

are 



(7 = at the moving boundary r = R , and C = Cq at f = oo , (4) 

and initial conditions are 

(7 = (7o and i? = i?o at t = , (5) 

where Cs denotes the solubility (or 'saturated mass concentration') of the 
solute in the fluid medium0and Co the initial uniform solute concentration. 

If we consider C{r,t) as a dimensionless variable ((7 — Cq)/{Cs — (7o), 
measure length in units of Rq and time in units of Rq/{ttD), the governing 
equations ([l])-(|5]) can be written in a nondimensionalized form 

dC ( p„\R^dRdC 1 d f odC\ , 



dt \ pmj dt dr nr"^ dr \ dr 

dR{t) (dC\ 



) , (7) 

r=R{t) 



dt \dr J 

C(i?(t),t) = 1 and C(oo,t) = 0, (8) 

C(r, 0) = and R{0) = 1 , (9) 
where t = TiDt/Rl, r = v/Rq, R = R/Rq, and e = {Cs — Co)/[npp{l — 



^Here the solubility Cs is treated as a constant, implying that the particle size effect 
on solubility as may be observed for submicron particles (often due to significant surface 
energy influence), is ignored for theoretical simplicity 



5 



Cs/ Pm)]- In genera equationsl (|6])-(9) describe a free-boundary (or moving- 
boundary) nonlinear problem, intractable to exact analytical solutions. But 
if e can be regarded as a small parameter (e.g., |e| << 1), which is generally 
valid for dilute solutions with materials of low solubility when Cs, C'o << Pp 
and Cs « Pm, the solutions to those nonlinear equations may be succes- 
sively approximated by solutions of linear equations following a perturbation 
procedure. 



3 Exact Quasi-Stationary Solutions 

Usually with perturbation approximation, the leading-order solution describes 
the most significant part of the phenomenon under study. Hence attention 
here is restricted only to the leading-order solutions. 

For small e, usually corresponding to the situation of relatively low sol- 
ubility, (7) indicates that the time variation of R{t) is slow comparing to 
that of C. Thus, the convection term in the convection-diffusion equation 
(6) can be neglected when considering the leading-order effect. In terms of 
perturbation solutions for small e, C can be written in the typical expansion 
form as 

C = C<°> + eC<'> + e^C<^> + ... , (10) 
where the zeroth-order solution C^°^ is the base solution at e = to the 
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zeroth-order equations 



dt 7rr^ dr \ dr J ~ 

C<°>(i?, t) = 1 and C<^>{oo, t) ^ , (13) 
C<o> (r, 0) = and i?(0) = 1 . (14) 



Obviously, at zeroth-order the particle radius R is not changing with time; 
it becomes a diffusion problem on a fixed domain R < r < oo. The solution 
of C^°^ to (11), which can also be written as 

d{rC<^>) _ 1 d\rC<°>) 
dt TT 

is given by 

C<«> = ^L-l^iy'^\-r,„l^ (15) 



r 



Noteworthy here is that treating as a time- independent "variable" in 
the zeroth-order diffusion equation is consistent with the so-called "quasi- 
stationary" approximation often used in engineering practice (where the con- 
vection transport is neglected and the diffusion equation is solved with the 
particle surface being considered stationary) . Here we ignore the initial con- 
dition R(0) — 1 and consider as a variable yet to be determined. 

To evaluate the dissolution process of a particle, we need to consider the 
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first-order effect of e at least for R{t) in (7) where only the zeroth-order C 



<o> 



given by (15) is involved, i.e., 

dR{t) (dC<'^>' 



dt 



dr 



1 1 

+ 



r=R{t) 



(16) 



As pointed out by Krieger et al. (1967), Chen and Wang (1989), and more 



recently [Rice and Do| ( |2006[ lp| ^ can be rearranged with some variable 
substitution to have a form of homogeneous ordinary differential equation as 



T- 



du 
dr 



e(2-e) + (w + e) 



or r 



u 



du 
dr 



1 + 



u 



(17) 



where r = u = R/r, i = y e/(2 — e), and u = {u + e)/^Je(2 — e) (which 
are mathematically valid for < e < 2). 

Straightforward integration of ( [l7| ) incorporating the initial condition 
R{0) = 1 yields 



(18) 



e(2-e)(l + n2) 

This is a solution that can only be expressed in an implicit form of R{t) with 



u 



{R/Vi + e)/ve(2 — e), though a cleaner formula than those presented 



■It seems though the analytical solution to (16) obtained by Krieger et al. (19671 was 



(1989 


) and 


Rice and Do 


(2006 


). However, 


Krieger et al. 



( 1967 ) used their "highly nonlinear" implicit formula merely to iteratively determine the 



value of diffusion coefficient, while Chen and Wang ( 1989 1 solved the same equation for an 



analytical solution with the application in drug particle dissolution in mind. The recent 



work of |Rice and Do ( 2006 ) again obtained "an exact analytical solution" by solving the 
same mathematical problem although with a slight change in the form of a parameter to 
account for the "bulk flow effect" as in (3)). 
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by previous authors (Krieger et al. 1967; Chen and Wang 1989 Rice and 



Do, 2006). The formula given by (18) suggests an easy way of generating 



curves for i? as a function of t by first selecting a series of values of m (> 



e) to calculate t (according to (18)), and then from the relationship R 



{uy/2 — e — y/e)\/~£t to calculate the corresponding R(t) from the given u and 



known t. Furthermore, as shown by previous authors (Chen and Wang, 1989 



Rice and Do, 2006), the time to complete dissolution to when i? = happens 



at M = e and is thus from (18) given by 



exp[e(2tan ^ e — tt)] exp[e(2tan ^ e — ti)] 



6(2-e)(l + e2) 



2e 



(19) 



Of mathematical interest, it would be worthwhile to mention that a seem- 
ingly different approach to the same mathematical problem was presented by 



Duda and Vrentas (1969) who were able to obtain the leading-order quasi- 



stationary solution of R{t) in an explicit form. With their approach, the 
moving interface is immobilized by introducing a boundary-fitted coordinate 
mapping 

r 

(20) 



r 



They also converted the equation system (6)-(9) into one in terms of a com- 
pound variable 



C = and C = C 



<o> 



eC 



<i> 



+ ... 



(21) 
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such that they arrived at 



dR{t) 
dt 



fR{t) \ df 



dlT 

Q(J<0> Q 

f 



71 dr'^ 



R{t) 



'dC 



1 < f < oo 



<0>~ 



df 



Similar to (15), the solution of (7*^°^ for (22) is then 



(22) 
(23) 



(:7<o>(l,t) = 1 and C'<°>(oo,t) = 0, (24) 
C'<o>(f, 0) = and i?(0) = 1 . (25) 



C<°> = 1 - ^ / e-'^'rfT^ and thus = -i . (26) 

V )f=i Vi 



Therefore, (23) yields 



R(t) = Jl-2e(2Vi + t 



(27) 



which is indeed a clean explicit formula. With this, Duda and Vrentas ( 1969 ) 



could also carry out derivations of subsequent higher-order perturbation so- 
lutions, which would be very difficulty, if not impossible, with the implicit 



formula (18). Now the question is why we can have two apparently different 



solutions (18) and (27) for the same order of approximation. To make sure 



(18) and (27) are reasonably equivalent mathematical solutions, it is helpful 



to check the numerical values with each of the formulas. Assuming e = 0.1 



(then e = 0.22942) and m = 1, (18) yields t = 1.83532 and then R = 0.45504 
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but ([27l) gives R = 0.30173. At m = e = 0.22942, 1^9^ predicts the time to 



complete dissolution to = 2.69711 when i? — )• 0, but the time to complete 



dissolution based on (27) would be [Jl + l/(2e) - if = 2.10102. Thus, we 



see that (27) cannot be the same as (18) at least for e = 0.1. 



A careful examination of the derivation of (22), however, reveals an error 
(which was somehow not corrected by those authors even in several follow-up 



publications, e.g., Duda and Vrentas, 1971 Vrentas et al. 1983). The correct 



expression of (22) should be 



d{C<^>) 1 d\C 



<0>^ 



dt 



1 < r < CO 



(28) 



Therefore, we should have 



C 



<o> 



2 l-ir-l)Ry/n/(it) 



TT JO 



e ^ drj and 



'dC 



<o>~ 



df 



f=i 



(29) 



which leads to the same equation as (16) and solution as (18) rather than 



(27). Thus, the same leading-order result can be obtained via seemingly 



different treatments. For describing the quasi-stationary dissolution process. 



(18) should be taken as the correct leading-order solution (for < e < 2). 

It might be noted that so far consideration is only given to the case 
of < e < 2, which describes the quasi- stationary dissolution process of a 
spherical particle. Mathematically, solution also exists for the case of e < as 



well as e > 2 in (16). From a physical point of view, the case of e < in (16) 



describes the inverse process of precipitation growth of a spherical particle. 



11 



i.e., when Cq > Cs corresponding to the situation of particle growth in a 



supersaturated solution. Somehow, the exact analytical solution to (16) for 



e < does not seem to have been presented in published literature, unlike 
the case of e > 0. Here it is derived to complete the mathematical solution 



for (16). With e < 0, (17) must be replaced by 



du [-e(2 - e)l - (u + e)^ du 

T— = or T— 

dr u dr 



u + i ' 



(30) 



where e = y — e/ (2 — e), and u = {u + e)/y—e{2 — e). The solution to (|30|) 
is then 



1 



[-e{2 - e)]{u^ - 1) \u-l 



u 



1 



(n > l,e < 0) 



(31) 



which appears to be quite different from (18)^ As with the dissolution 
case, curves of R{t) can easily be generated by first selecting a series of 



values of u (> 1) to calculate t from (31), and then from the relation R 
e + y^^)y/—et to calculate R(t), for the particle growth case. 



{uV2 



Although the quasi-stationary diffusion-controlled mass transfer problem 
can be treated as a leading-order perturbation problem, it is noteworthy that 
the perturbation procedure is only performed on the convection-diffusion 



■^Actually if the identity ln[(u + l)/{u — 1)] = 2 coth is used, we can have (31 1 in a 
similarly looking form to ( 18 1 as 



2 ^ exp (2e coth ^u) 



({t > l,e < 0) 



-e(2-e)](u2-l) 
where coth~^ denotes the inverse hyperbolic cotangent function. 
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equation (6) in terms of the small parameter e such that the neglection of 
the convection term can be justified in the zeroth-order equation. In fact, 
the convection term may naturally disappear in (6) when pp/ — )■ 1 (which 
might in fact be a quite reasonable assumption especially for a solid solute 
particle either dissolving in a liquid solvent or growing by precipitation in a 
supersaturated liquid solution, because the density of solid solute is usually 
not much different from that of liquid solution, unlike the situation of a liquid 
droplet in gas or a gas bubble in liquid with orders of magnitude of density 
differences). In that case, the perturbation treatment in terms of e becomes 
unnecessary; the "leading-order" solution is the exact solution for any value 



of e. Then, (30) and (31) with e being replaced by — e are also valid for the 



case of e > 2 but with R = {uy/e — 2 — ^/e)^/ei which becomes zero when 



u = i, where t ^ \J e/ {e — 2) (> 1) and u = [u + e) / \J e{e — 2)^ for u> e. 
Special attention though should be paid at e = 2 where 

du {u + 2Y 

^-r = ' y^^) 

ClT u 
which has the (implicit) solution 



exp[-2/(^ + 2)] exp[-A/iu + 2)] ,^ . ^ . 

^2 °^'= iu + 2Y (0<«<oo), (33) 



with R = UT = u\/i. As expected, (18) indeed approaches (33) at the limit 



as e — 7- 2 by applying L'Hospital's rule, and so does (31 ) with e being replaced 
by — e for e > 2 by using the relationship e = limj;_j.o(l + x)^/^. 
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Thus, we have a complete set of solutions to (16) with (18) for < e < 



2, (31) for e < and e > 2 (with e being replaced by — e), and (33) for 
e = 2, covering all possible e 7^ (whenever ignoring the convection term 
in (6) is justifiable in practice). Moreover, the case of e = corresponds 
to a mathematically trivial solution = 1, if desired to be included for 
completeness. 



4 Approximate Explicit Formulas 



Due to the awkwardness of practical usage of the implicit solution to (16) 



authors (e.g., Chen and Wang 1989 Rice and Do 2006) who derived the 



analytical solution (for < e < 2) often also suggested explicit "approximate 
solution" with the quasi-steady-state treatment (when the time derivative 
term in the diffusion equation (11) is also neglected). Simple as it may 
appear, however, the explicit quasi-steady-state solution may be found not 



to offer a satisfactory approximation to (18) for e > 10~ (depending on 
the intended application). To provide much improved approximate explicit 
formulas for practically accurate evaluation of R{t), effort is made here via 
a semi-empirical approach. 

In view of the physical meanings, the two terms on the right side of 



(16) represent two different aspects of the diffusion process: 1/R comes 
from the spherically symmetric steady-state solution of the Laplace equation 
whereas describes transient diffusion from a planar surface. Initially 
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when \/i « R, the concentration gradient front has not propagated far 
enough from the particle surface to bring the curvature effect out yet and 
thus the effect of transient diffusion from a planar surface dominates. With 
increasing t, at some point \/i » R is expected to happen especially for 
|e| << 1; then the steady-state diffusion term becomes dominant and (16) 
can be reduced to 

dR(t) _ . 



dt R{t) ' 

which is a so-called "quasi-steady-state" approximation commonly seen in 



Duda and Vrentas , 


1971 


Vrentas et al. 


1983 


Chen and 



Wang[ 1989; Rice and Do, 2006). The quasi-steady-state solution to (34) is 

R{t) = Vl-2et foit»R\ (35) 



simply 



which yields noticeable difference from (18) even at e = 10 ^. 
On the other hand, for t << R^ we have 



dR{t) 



and thus R{t) = 1 - 2eVi . 



(36) 



dt y/i 

Intuitively then, one might believe that a constructed explicit formula like 

R{t) = Vl - 2et - 2eVi (37) 
could provide an improved approximation from the quasi-steady-state result 



(35) to the exact solution (18). Clearly, (37) approaches R = 1 — 2e\/i as 
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Table 1: Comparison of predicted values of time to complete dissolution from 
the quasi-stationary, quasi-steady-state, and intuitively constructed models, 
with relative error with respect to to given in the parentheses. 



to from (19) l/(2e) from psf to from (38) 



1 


0.1039 


0.5 (381%)~ 


0.1667 (60:4%) 


0.5 


0.2984 


1 (235%) 


0.5 (67.6%) 


0.1 


2.6971 


5 (85.4%) 


4.1667 (54.5%) 


0.05 


6.3622 


10 (57.2%) 


9.0909 (42.9%) 


0.01 


40.421 


50 (23.7%) 


49.020 (21.3%) 


0.005 


85.876 


100 (16.4%) 


99.010 (15.3%) 


0.001 


466.54 


500 (7.2%) 


499.00 (7.0%) 


0.0005 


952.01 


1000 (5.0%) 


999.00 (4.9%) 


0.0001 


4890.6 


5000 (2.24%) 


4999.0 (2.22%) 



t — )■ (where the 'higher-order' term associated with t = (y/i)'^ becomes 
negligible). The approximate time to complete dissolution corresponding to 



(37) is 



1 



(38) 



°"2e(2e+l)- 

For comparison at various values of e, the values of time to complete dis- 
solution are tabulated in table 1 for to from the exact (quasi- stationary) 



solution (19), the quasi-steady-state result l/(2e), and the intuitively con- 



structed model result to from (38). The consistent improvement of (38) over 



l/(2e) from (35) in comparison with (19) is obvious, but not as significant 



as desired. Both (35) and (37) underestimate the change of R with t. 



Interestingly, (27) obtained by Duda and Vrentas (1969) is somehow found 



to offer a better approximation than (37) to the exact quasi- stationary solu- 
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tion especially for e < 0.1, although the its derivation has been shown not 
quite correct in the mathematical sense. For example, the time to complete 



dissolution predicted with [Jl + l/(2e) — 1]^ based on (27) would be 0.0505 



2.1010, 37.717, and 457.23 respectively for e = 1, 0.1, 0.01, and 0.001. More- 



over, (27) tends to predict faster dissolution whereas (37) slower. This is 



because (27) over estimated the flux term associated with 1/yt in (26) by 



mistakenly replacing R (< 1) with 1. However, the effect of l/y/i usually 
only dominates for a short time when t is small and R is not too far from 
unity especially when e << 1. Shown in Fig. 1 is a comparison among the 
exact quasi-stationary solution (18), the quasi-steady-state solution (35), and 



the approximate formulas (27) of Duda and Vrentas (1969) and (37). Even 



at e = 0.01, the deviation of quasi-steady-state solution (35) from the exact 
solution is still quite significant due to the unaccounted initial effect from 
the flux term l/\/t for small t. The overall improvement of (37) from the 
quasi-steady-state solution (35) is clear, especially for small t (or in general 
t < 0.5to) where the curve of (37) consistenly remains close to that of the 
exact quasi-stationary solution. 

In view of Fig. 1, an explicit approximation formula may be constructed 



semi-empirically by combining (27) and (37) as 



ale) 1 - 2e(2v^ + t)] + [1 - a(e)] (Vl -2tt- 2eVi 



(39) 
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Figure 1: Comparison among the exact quasi-stationary solution (18), the 



quasi-steady-state solution (35), and the approximate formulas (27) of Duda 
and Vrentas (1969) and (37) for dissolution of particle at e = 0.01. 
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which with 



a e 



0.781 [1 - 1.935/ (1 + 1.05 g-o-^^TS)] (o < e < 0.1) 
0.0193(logio tf - 0.2703 logio e + 0.095 (0.1 < e < 0.5) 



can consistently produce the value of R{t) very close to that of the exact 



solution (18) for any e < 0.5 in the entire range of < t < to- By the similar 



token, accurate formulas of a{e) for e > 0.5 can also be constructed but is 
not attempted here because most practical situations, e.g., in pharmaceutical 



dosage form development (Curatolo, 1998 Kerns and Di 2008), typically 
concern with e < 0.1. 

For particle growth in a supersaturated environment, i.e., Cq > Cs as in 
the case of phase separation in a drying coating, the exact quasi-stationary 



solution is given by (31) with e < 0. Fig. 2 for particle growth shows that 



comparing to (31), (37) seems to overestimate the particle growth whereas 



(27) of Duda and Vrentas (1969) and the quasi-steady-state solution (35) 



both underestimate it. 

Based on this observation, a fairly accurate explicit approximate for- 



mula to the exact quasi-stationary solution (31) may be semi-empirically 



constructed with the same form as (39) but having 



a e 



0.5381 



0.3 



1 + |e| 



-0.6514 



;-0.5 < e < 0) , 



(40) 



for evaluating the situation of particle growth. 
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Figure 2: Comparison among the exact quasi-stationary solution (31), the 



quasi-steady-state solution (35), and the approximate formulas (27) of Duda 



and Vrentas (1969) and (37) for particle growth by precipitation at e = 0.01. 
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5 Concluding Remarks 

Starting from a consolidated mathematical formulation of the spherically 
symmetric mass-transfer problem, the quasi-stationary approximating equa- 
tions can be derived based on a perturbation procedure for the leading-order 
effect. For the diffusion-controlled quasi-stationary process, a mathemati- 
cally complete set of the exact analytical solutions is obtained in implicit 
forms for consideration of both dissolution of a particle in a solvent and 
growth of it by precipitation in a supersaturated environment. Understand- 
ing the dissolution behavior of solid particles in liquid plays an important 



role in pharmaceutical dosage form development (Chen and Wang, 1989 



Rice and Do 2006), and particle growth by precipitation in a supersaturated 



environment is relevant to the drug-polymer microsphere formation process 



(Wu, 1995) as well as the observed phase separation process during solvent 



removal in drying of a coating with the drug-polymer mixture (Barocas et 



al. 2009 Richard et al. 2009). The commonly used explicit formula based 
on the solution with quasi-steady-state approximation is shown to provide 
unsatisfactory accuracy unless the solubility is restricted to very small values 
(corresponding to e < 10~^). Therefore, accurate explicit formulas for the 
particle radius as a function of time are also constructed semi-empirically to 
extend the applicable range at least to —0.5 < e < 0.5 for practical conve- 
nience. 
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